This files shows the basic steps of fly eye measurement and feature extraction with one image as an example.

# Load packages 
library(lattice)
library(ggplot2)
library(sp) 
library(raster) 
library(tiff)
library(EBImage)
library(Gmedian)
library(gridExtra)


# source helpful R code and specify image folder path 
source('./src_Fly.eye.pat_funcs.R')
dir.in='./example_images/'
dir.out='./output/'
cat("check dir: ",dir.in,"; dirout",dir.out,"\n");
## check dir:  ./example_images/ ; dirout ./output/
# Read in images and begin process
image.filenames <- list.files(path=dir.in,pattern="*jpg$", full.name=F)
image.filenames=image.filenames[1]; #use one image as an example
print(image.filenames);
## [1] "test1.jpg"
n.images=length(image.filenames);
image.files=paste(dir.in,image.filenames,sep='/');

# process one image at a time
# save each image quantitive features in a lsit

  image.i = 1
  image = readImage(image.files[[image.i]] )
  image.name = image.filenames[image.i]
  cat('begin processing image',image.name,',',image.i,'images out of',n.images,'total images\n')
## begin processing image test1.jpg , 1 images out of 1 total images
  display(image, method="raster")

  dim(image)
## [1] 2560 1920    3
  # Resize to fit memory
  image.res=resFunc(image)
  display(image.res, method="raster")
  dim(image.res)
## [1] 640 480   3
  # Assign resized images RGB channels to data frames
  image.Ori <- RGBintoDF(image.res)
  head(image.Ori)
##   x   y R G B
## 1 1 480 1 1 1
## 2 2 480 1 1 1
## 3 3 480 1 1 1
## 4 4 480 1 1 1
## 5 5 480 1 1 1
## 6 6 480 1 1 1
  # White TopHat morphological transform
  image.Top <- wTopHat(image.res,y=5,z='diamond')
  
  # Display images
  par(mfcol=c(1,2))
  dispImg(image.res)
  dispImgT(image.Top,0.99)

  ## Threshold and centroids
  # Assign gray channel to data frame
  image.grey <- GintoDF(image.Top)
  head(image.grey)
##   x   y G
## 1 1 480 0
## 2 2 480 0
## 3 3 480 0
## 4 4 480 0
## 5 5 480 0
## 6 6 480 0
  # Threshold to keep pixels with intensity > 0.99 quantile
  image.thres.retain<-image.grey[image.grey$G > quantile(image.grey$G,0.99), ]
  
  # select ROI, region of interest through two rounds of pixel fitering based on distances to centroids
  image.thres = image.thres.retain
  cutoffs=c(0.8,0.3);
  for(cutoff in cutoffs){
    cat('ROI selection with cutoff',cutoff,'\n')
    # Estimate images centroid
    centroids <- Weiszfeld(image.thres[,1:2])
    
    # extract each pixel x,y coordiantes retained pixels
    thresXY <- image.thres[,1:2,drop=FALSE]
    p1<-ggplot(thresXY,aes(x=x,y=y))+geom_point()+
      geom_point(aes(x=centroids$median[1], y=centroids$median[2]),colour="red",size=4)+
      theme_bw()+ ggtitle(image.name)+xlab('')+ylab('')+
      theme(panel.grid = element_blank(),
            axis.text = element_blank(),
            axis.ticks = element_blank())
    ## Distances to centroide,calculate distances to centroid
    pdist <- pointDistance(p1=thresXY, p2=centroids$median,lonlat=F)
    
    # Mark distances > 0.8 quantile, as they belong to
    # points outside the eye boundary in their majority
    pLogic <- (pdist < quantile(pdist, cutoff));
    distCent <- cbind(distCent = pdist, selected = pLogic)
    
    # Plot example histograms
    p2<-ggplot(data.frame(distCent),aes(x=distCent))+geom_histogram(bins = 30)+
      geom_vline(xintercept = quantile(distCent[,1], cutoff),col="red",lty="dashed", lwd=2)+
      theme_bw()+ ggtitle(image.name) 
           
    
    # Join thresholded and distances lists
    thresDist <- cbind(image.thres,distCent);
    # retain points with distance < quantile cutoff 0.8
    distSelect <- thresDist[thresDist$selected==1,]
    
    # Plot examples (black clouds with blue roi overlay)
    p3 <- ggplot(data=thresDist, aes(x=x, y=y,color=selected))+
        geom_point(show.legend = FALSE) + 
        theme_bw()+ ggtitle(image.name)+xlab('')+ylab('')+
      theme(panel.grid = element_blank(),
            axis.text = element_blank(),
            axis.ticks = element_blank())
   
    print(grid.arrange(p1,p2,p3,ncol=3))
   
    # Overlay to image and plot example
    p <- ggplot(data = image.Ori, aes(x = x, y = y)) +
        geom_point(colour = rgb(image.Ori[c("R","G","B")])) +
        labs(title = "Original Eye selected Points",
             cex=0.5) +xlab("x") +ylab("y") +
        geom_point(data=distSelect, alpha=0.2) +
        plotTheme()
    
    image.thres <- distSelect[,1:3]
  }
## ROI selection with cutoff 0.8

## TableGrob (1 x 3) "arrange": 3 grobs
##   z     cells    name           grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]
## 3 3 (1-1,3-3) arrange gtable[layout]
## ROI selection with cutoff 0.3

## TableGrob (1 x 3) "arrange": 3 grobs
##   z     cells    name           grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]
## 3 3 (1-1,3-3) arrange gtable[layout]
  ## Subset and confidence ellipse and Add ellipse to plot 
  ellPlots = ellPlot(distSelect,0.90)
  ellPlots #it's a ggplot object

  # Extract components
  build = ggplot_build(ellPlots)$data
  ell = build[[2]]
  
  # Select original image points inside ellipse
  orig.ell = data.frame(image.Ori[,1:5],
             in.ell = as.logical(point.in.polygon(image.Ori$x,image.Ori$y, ell$x, ell$y)))
  orig.pix = orig.ell[orig.ell$in.ell==TRUE,]
  
  # plot example: raw eye image + ROI circled in blue
  p <- ggplot(data = image.Ori, aes(x = x, y = y)) +
    geom_point(colour = rgb(image.Ori[c("R", "G","B")])) +
    labs(title = "Original Eye final ROI", cex=0.5) +
    xlab("x") + ylab("y") +
    geom_polygon(data=ell[,1:2], alpha=0.2,
                 size=1, color="blue") +#plotTheme()
    theme_bw()+xlab('')+ylab('')+
      theme(panel.grid = element_blank(),
            axis.text = element_blank(),
            axis.ticks = element_blank())
  print(p)

  ## Create image from final ROI
  rois.image = roitoImg(orig.pix)
  display(rois.image,method='raster'); #check image
  
  ## begin extract ommatidia configuration features
  pic=rois.image
  #hist(pic);grid()
    
  y = equalize(pic)
  #hist(y);grid()
  #display(y, title='Equalized Grayscale Image',method='raster')
    
  grayimage<-channel(y,"grey")
  #display(grayimage,method='raster')
  nmask = thresh(grayimage, w=5, h=5, offset=0.02); #display(nmask)
  nmask = opening(nmask, makeBrush(3, shape='disc')); #display(nmask)
  nmask = fillHull(nmask); #display(nmask)
  nmask = bwlabel(nmask)
  cat("Number of ommatidia,",max(nmask),"\n");
## Number of ommatidia, 75
  fts = data.frame(computeFeatures.moment(nmask))
  head(fts)
##       m.cx      m.cy m.majoraxis m.eccentricity    m.theta
## 1 50.50000  1.500000    4.472136      0.8944272  0.0000000
## 2 38.06452  5.903226    8.830081      0.8437787  1.5130830
## 3 62.42623  5.459016   18.670999      0.9729211  0.3624306
## 4 45.89474  6.421053    5.544698      0.6322809 -0.2237600
## 5 31.25926  9.925926    6.518058      0.5827875 -0.1450000
## 6 54.40741 10.925926    6.421501      0.5724126  0.2180192
  fts2 <- data.frame(computeFeatures.shape(nmask))
  head(fts2)
##   s.area s.perimeter s.radius.mean s.radius.sd s.radius.min s.radius.max
## 1      8           8      1.144123   0.4370160    0.7071068     1.581139
## 2     31          19      2.759500   0.8339106    1.1263650     3.846429
## 3     61          35      4.752899   2.2218203    1.5452363     8.321352
## 4     19          13      1.986977   0.3333103    1.3931362     2.560510
## 5     27          16      2.465936   0.3606885    1.9465514     3.006503
## 6     27          15      2.471452   0.3104342    1.8961950     2.988125
  par(mfrow=c(1,5));

  display(pic,"raster")
  display(y, title='Equalized Grayscale Image',method='raster')
  display(grayimage,method='raster')
  display(nmask,"raster");
  display(nmask,"raster");
  text(fts[,"m.cx"], fts[,"m.cy"], 
       labels=seq_len(nrow(fts)), col="red", cex=0.8)

  par(mfrow=c(1,2));
  display(pic,"raster")
  #display(y, title='Equalized Grayscale Image',method='raster')
  #display(grayimage,method='raster')
  display(nmask,"raster");

  #display(nmask,"raster");
  #text(fts[,"m.cx"], fts[,"m.cy"],labels=seq_len(nrow(fts)), col="red", cex=0.8)
    
  # fts and fts2 contain raw measurements on each ommatidium, which can be written out to a file
  #out.file1=paste(dir.out,"/",image.name,"-coord.txt",sep='');
  #out.file2=paste(dir.out,"/",image.name,"-area.txt",sep='');
  #write.table(fts,file=out.file1,quote=F,row.names = F);
  #write.table(fts2,file=out.file2,quote=F,row.names = F);
  
  ## ROI seleciton done, begin to compute features based on raw ommadium measurements
  ## some filter:
  # images contain less then 6 ommaditia, discard
  # when calculating each feature, remove some extreme values, the largest 3 and smallest 3.
  # number of ommatida
  n.ommatidia=nrow(fts)
  if(n.ommatidia<6){ result[[image.name]]=NA;next} 
  # There are four features computed from data.frame fts 
  # nn: pairwise ommatidia nearest neighbor distances
  # cc: each ommatidium essentricity
  # nn.mean, nn.sd, cc.mean, cc.sd
  colnames(fts)[c(1,2,4)]; #x,y coordinates, essentricity
## [1] "m.cx"           "m.cy"           "m.eccentricity"
  nn.dist=c();
  for(i in 1:nrow(fts)){
    pairwise.dist=c();
    for(j in 1:nrow(fts)){
      if(i==j){next}
      distance=((fts[i,]$m.cx-fts[j,]$m.cx)^2+(fts[i,]$m.cy-fts[j,]$m.cy)^2)^0.05;
      pairwise.dist=c(pairwise.dist,distance);
    }
    nn.dist=c(nn.dist,min(pairwise.dist))
  }
  #length(nn.dist)==nrow(fts)
  nn.out=get_mean_sd(nn.dist)
  names(nn.out)=c('nn.mean','nn.sd')
  
  cc.out=get_mean_sd(fts$m.eccentricity)
  names(cc.out)=c('cc.mean','cc.sd');
  
  x0=c(n.ommatidia,nn.out,cc.out)
  names(x0)=c('n.ommatidia','nn.mean','nn.sd','cc.mean','cc.sd')
  
  # There are 12 features computed from data.frame fts2 
  # mean and sd for each column variable
  colnames(fts2)
## [1] "s.area"        "s.perimeter"   "s.radius.mean" "s.radius.sd"  
## [5] "s.radius.min"  "s.radius.max"
  out=apply(fts2,2,get_mean_sd)
  x1=out[1,]
  x2=out[2,]
  name1=paste0('mean(',names(x1),')')
  name2=paste0('sd(',names(x1),')')
  names(x1)=name1;names(x2)=name2
  out=c(x0,x1,x2)
  
  plot.new()
  grid.table(data.frame(out))

sessionInfo()
## R version 3.6.1 (2019-07-05)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS High Sierra 10.13.6
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] gridExtra_2.3   Gmedian_1.2.4   EBImage_4.28.0  tiff_0.1-5     
## [5] raster_3.0-7    sp_1.3-2        ggplot2_3.3.2   lattice_0.20-38
## 
## loaded via a namespace (and not attached):
##  [1] tidyselect_1.1.0    locfit_1.5-9.1      xfun_0.11          
##  [4] purrr_0.3.3         colorspace_1.4-1    vctrs_0.3.2        
##  [7] generics_0.0.2      htmltools_0.4.0     yaml_2.2.0         
## [10] rlang_0.4.7         pillar_1.4.6        glue_1.4.2         
## [13] withr_2.1.2         BiocGenerics_0.32.0 jpeg_0.1-8.1       
## [16] lifecycle_0.2.0     robustbase_0.93-5   stringr_1.4.0      
## [19] munsell_0.5.0       gtable_0.3.0        htmlwidgets_1.5.1  
## [22] codetools_0.2-16    evaluate_0.14       labeling_0.3       
## [25] knitr_1.26          parallel_3.6.1      DEoptimR_1.0-8     
## [28] Rcpp_1.0.3          scales_1.1.1        abind_1.4-5        
## [31] farver_2.0.3        RSpectra_0.15-0     png_0.1-7          
## [34] digest_0.6.22       stringi_1.4.3       dplyr_1.0.2        
## [37] grid_3.6.1          tools_3.6.1         bitops_1.0-6       
## [40] magrittr_1.5        RCurl_1.95-4.12     tibble_3.0.3       
## [43] crayon_1.3.4        pkgconfig_2.0.3     ellipsis_0.3.0     
## [46] MASS_7.3-51.4       Matrix_1.2-17       rmarkdown_1.17     
## [49] R6_2.4.1            fftwtools_0.9-8     compiler_3.6.1
#installed.packages()[names(sessionInfo()$otherPkgs), "Version"]